Self-organized critical network dynamics 
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We propose a simple model that aims at describing, in a stylized manner, how local breakdowns 
due unbalances or congestion propagate in real dynamical networks. The model converges to a self- 
organized critical stationary state in which the network shapes itself as a consequence of avalanches 
of rewiring processes. Depending on the model's specification, we obtain either single scale or scale- 
free networks. We characterize in detail the relation between the statistical properties of the network 
and the nature of the critical state, by computing the critical exponents. The model also displays a 
non-trivial, sudden, collapse to a complete network. 

PACS numbers: : 05.40.-a, 45.70.Ht, 89.75.Fb 
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Complex networks underlying many social and tech- 
nological systems is a subject of booming recent interest 
P, U • On one side the structure of such networks has 
non-trivial properties 0, which influences dramatically 
the nature of processes taking place on them (see e.g. 
ISb^J)- On the other, the network's structure constrains 
in a pecuhar way the growth and evolution 0, Q of 
the network itself. This calls for an extension of statisti- 
cal physics, which traditionally studies collections of dy- 
namical variables interacting through a fixed network, to 
systems where the network of interactions itself becomes 
a dynamical variable. 

Here we focus on dynamical networks where links do 
not represent physical bonds but rather relationships or 
communication channels. Ref. j3| suggests that a struc- 
tured network of communications aimed at solving prob- 
lems or carrying out specific functions is a crucial feature 
of firms and organizations in general. Routers table in the 
Internet is also an example of a dynamic communication 
network. The network of financial institutions, linked by 
mutual contracts and loans, provides a further example 
of a dynamic network. Beyond "static" design problems, 
such as e.g. minimizing congestion 3] or redistributing 
optimally the loads 8], systems of this type also pose 
"dynamical" problems such as how and to what extent 
does congestion or breakdown events propagate through 
the system. 

Here we address these problems in a dynamic network 
subject to two competing forces: on one side, there is a 
drive toward increasing complexity, by e.g. adding new 
links, because the system performs more efficiently its 
functions as it becomes more densely interconnected. On 
the other, the resulting increase in complexity may bring 
about conflicting constraints, imbalances or congestion 
problems which may cause a local breakdown of the net- 
work. A local breakdown may engender a re-adaptation 
in its neighborhood which may inadvertently cause the 
breakdown to propagate further on the network. 

For example, a change in some router's table in the In- 
ternet, which is meant to improve efficiency, may inadver- 
tently cause congestion at some node downstream. This 
may trigger several other changes in that local neighbor- 
hood, as routers try to avoid the congested nodes. But 



these changes may, in their turn, cause further conges- 
tion elsewhere, and the problem may expand even fur- 
ther, as an avalanche, to a wider region. The stipulation 
of a contract between two institutions, which is in prin- 
ciple beneficial to both, may also increases their opera- 
tive constraints, making them less adaptable to a chang- 
ing environment and hence more exposed to the risk of 
bankruptcy. The failure of one institution, likely induces 
a rearrangement of the institutions linked with it and 
perhaps engenders effects which propagate further across 
the network [^. Similar phenomena may take place in 
social or trade networks. 

Rather than trying to model in a realistic manner one 
of the problems just discussed, we focus on a simple 
model of network dynamics which captures the two main 
ingredients discussed above: a slow dynamics where links 
are added to the network and a fast relaxation dynam- 
ics of avalanche events. The motivation for this choice, is 
that, the detailed understanding of the behavior of a sim- 
ple model with these features, may be the basis or at least 
a guide for addressing more complex and realistic situa- 
tions, such as those discussed above. Our main finding 
is that such systems can self-organize close to a critical 
point where each modification of the network's architec- 
ture can have unforeseable consequences which possibly 
affect a wide region of the system. This may have some 
bearing on the intermittence of Internet traffic or on the 
nature of financial crises and recessions. 

We start from a empty network of N nodes in which 
every node i is assigned a fitness fi drawn from a proba- 
bility distribution p{f ). Let Fi be the set of neighbors of 
i and ki = \Fi \ be the number of neighbors oii. At every 
time step a link is added between two previously uncon- 
nected random nodes i and j ^ Fi. With probability fj 
nothing happens whereas with probability 1 — fj the node 
j becomes unstable or congested and it "topples" . As a 
result all its links (including that with i) are rewired to 
randomly chosen nodes, i.e. for any h € Fj, a node I F^ 
is chosen at random and the link jh is rewired to hi. In its 
turn, with probability 1 — fi, also node I may become un- 
stable and topple. Hence toppling of node j may start an 
avalanches of toppling events which propagates through 
the network rearranging it. Unstable nodes, after they 
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topple, remain unconnected from the network and 
are assigned a new fitness value drawn from the distribu- 
tion p{f). Hence toppling is equivalent to replacing the 
unstable node with a new one. 

In order to stabilize the network and reach a stationary 
state, we introduce dissipation of links: at each toppling 
event, with probability A all the links of the unstable 
node are removed from the network. Note that without 
dissipation the number of links would increase in time 
until the complete graph is reached. The complete graph 
(j G Fi i ^ j) is an absorbing state of the dynamics, 
because no link can be added to it. 

The distribution p{f) is the only parameter of the 
model. An alternative class of models can be defined by 
specifying the probability that a node with k neigh- 
bors becomes unstable upon addition of a further link. 
A relation between the two models is possible, along the 
lines of Ref. the limit A —> on no dissipation. 

Then the probability to find a node with k neighbors and 
h e [/,/ + df) is p{f\k)df^^ fpif)df where /'= is the 
probability that a node with fi = f has k neighbors. 
Then a model with 
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is completely equivalent to one specified in terms of p{f), 
in the limit A ^ 0. The dependence on k of Uk reflects 
the fact that ki and fi are positively correlated because 
nodes with higher fitness have smaller chance of becom- 
ing unstable. For convenience we shall refer mostly to 
models specified in terms of Uk using Eq. to translate 
the results in the original model. 

Our results can be summarized as follows: i) when Uk 
decays faster than 1/fc there is a critical Ac such that 
the network evolves toward a complete graph for A < Ac. 
The same happens for Uk — b/k {k ^ 1) and b < 3/2 
whereas when b > 3/2 or when Uk decays slower than 
1/fc, the collapse takes place only in finite networks. In- 
deed we find complete graphs only for A < Ac ^ N~'', 
where 7 — ^^5r for 3/2 < 6 < 2 and 7=1 otherwise. 
ii) the non-collapsed phase A > Ac is characterized by an 
uncorrelated random network |12| with finite average de- 
gree and a degree distribution pk which depends on p{f) 
(or Uk) (see Fig. In particular, if Uk decays slower 
than 1/fc then Pk decays faster than any power, whereas if 
Uk ^ b/k then pk ~ k^^. Hi) The dynamics converges to 
a stationary sequence of avalanches of rewiring processes 
with a power law distribution P{s) s~'^ of sizes (see 
Fig. EJ. As in Ref. [ij, we shall define the size s of an 
avalanche as the number of toppling events that it causes. 
The exponent r takes the mean field value[l3 r = 3/2 
when Uk decays slower than 2/fc, whereas r = 1 + 1/6 if 
Uk ^ b/k with 3/2 < 6 < 2. 

The collapse to a complete graph, where congestion is 
minimal, is reasonable, given that the network is trying 
to adapt by avoiding congestion. The non-trivial issue 



on which we shall concentrate mostly, is related to the 
self-organized critical state. In view of the special role 
played by the case Uk ~ 1/fc, we focus on the following 
simple forms of p{f) and to the corresponding Uk'. 
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In order to shed light on the model's behavior, let us 
derive an equation for the avalanche distribution. Define 
Sk as the avalanche size originating from an unstable node 
with fc neighbors. This is a random variable which can 
be decomposed as follows 
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Sfc = 1 + d ^ Vkj Skj ■ 
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into the contribution of the unstable node and those of 
the avalanches Sk^ ensuing from its neighbors, with kj 
being the number of neighbors of the neighbor. Note 
that the sum runs over fc + 1 links as it includes the 
link which caused the instability and the fc pre-existing 
neighbors. In Eq. © d = 0, 1 describes the effect of 
dissipation with P{d = 0) = A = 1 — P{d = 1) whereas 
Vkj takes value Vk^ = if rewiring of the link to the 
j'^ neighbor causes no further toppling, and Vkj = 1 
otherwise. Hence P{vi = 1) = 1 — P{vi = 0) = ui. 

Now we can write the generating function (pkiz) ~ 
{z^'') of the probability P{s\k) to have an avalanche of 
size s given that the initiator node has fc neighbors. From 
this, it is easy to find the generating function xi^) of fhs 
distribution P{s) of avalanche sizes 

-. 00 

x(z)^^P(s)z^ = ^^p,.z.,</>fc(z). (4) 

S fc = 

with u = '^f.pkUk- After some algebra, using the fact 
that the rewired nodes are chosen randomly in the net- 
work, Eq. Q leads us to 
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which is a non-linear self-consistent equation for xi'A- 
Note that x(l) = 1 as it should and that 
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x'(i) = = Y^TTT 



(l-A)((fc + lW 



Apart from the parameters A and Uk, Eq. |(SJ de- 
pends also on the stationary state degree distribution pk ■ 
Hence, before discussing Eq. further, we need to elab- 
orate on the nature of the stationary state. A necessary 
condition in order to be at the stationary state is that 
the total number of links K(t) present in the network is 
constant on average. K{t) increases by one for the ran- 
dom addition of a new link and is reduced by the total 
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FIG. 1: Degree distribution pk for network of A'^ = 10* nodes 
(averaged over SOOA'^ time steps after equilibration) for the 
model of Eq. Q and 6 = 1.8, . . . , 3.4 and A = 5 • 10"^ and 
for Uk = u — 0.5 (inset). The solid lines corresponds to the 
theoretically expected pk in the limit A ^ [Eq. 1101 and 
Pk = ^(1 — u)'^ respectively]. 

number of links A dissipated during the avalanche that 
follows. Here 

A^vJ2d.{h + l) (7) 

1=1 

where w = 1 if the chosen site is unstable and v — 
otherwise, and di — 1 if dissipation occurs at the toppling 
site i, otherwise di = 0. Consequently v and di have 
average values u and A and thus (A) = A ((fc + l)uk) (s). 
Then, stationarity (A) = 1 and Eq. © imply that 

{{k + l)uk)^l. (8) 

As a byproduct, we find (s) = 1/A, in perfect agreement 
with numerical simulations. 

The stationary degree distribution pk can be found by 
quantifying the Markov chain of possible transitions ki — > 
kl during the dynamics. In the limit A — > and N —^ oo 
where we can neglect dissipation and finite size effects, 
there are only two processes which take place on each 
node: fc^ — » fc^ + 1, with probability 1 — Wfc. and ki —> 
with probability Uki ■ Then, in the stationary state, pk 
satisfies 

Pk+i = (1 - Uk)pk- (9) 

Taking the sum of Eq. on k we find po — (ufc), which 
means that the fraction of sites with no neighbor is equal 
to the probability that a node becomes unstable. Fur- 
thermore multiplying Eq. by fc + 1 and taking the 
sum over fc, we recover the stationary condition ||SJ). In 
the simplest case Uk = u for all k, we find pk — u{l — u)^ 
whereas with Eq. (j^J we find 

where the asymptotic power law behavior holds for k ^ 
1. Notice that (fc) = !/(& — 2) diverges when & — > 2~ and 
that there is a finite fraction u = 1 — 1/6 of unconnected 
nodes. Still the network has a giant connected compo- 
nent for b < 7 [H. Eq. © yields a Pk which decays 
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FIG. 2: Average degree as a function of A for different system 
sizes and b = 1.8. The value of Ac at which the transition 
takes place is plotted in the inset against TV. The full line is 
the theoretical prediction Ac oc N''' with 7 = 0.75 for b = 1.8. 



faster than any power if Uk decays less slowly than 1/fc, 
or if it increases. Conversely, ifuk decays faster than 1/fc, 
we find that pk is not normalizablc for N 00. Numer- 
ical simulations (see Fig. ^| fully supports this picture, 
even though the effects of dissipation and finite size are 
clearly evident for fc ^ 1. 

The neglect of dissipation, when N is finite, is a rea- 
sonable approximation if nodes with maximal degree 
ki = N — \ are not stable. A node connected to all 
neighbors, cannot receive further links and hence cannot 
become unstable. Its degree decreases only if dissipation 
occurs at a node connected to them. The rate of this 
process, for a node with degree fc, is Afc/ (fc). Hence if 
\N / (fc) ^ 1, nodes with k ^ N decay very fast and the 
only effect of dissipation is to introduce a cutoff fee ^ 1/A 
in the distribution pk- When \N/ (k) ~ 1 we expect a 
transition (close) to the complete graph ki — N — 1 for 
all i. When {k) is finite, i.e. for 6 > 2 or for Uk which 
decays slower than 1/fc, the collapse to a complete graph 
takes place for A < Ac N-^. When 3/2 < 6 < 2 

2-b 

the average degree [k) ^ Ni>-^ diverges with the system 
size. Then the decay rate of totally connected nodes is 
^ XN'' with 7 = ^'^d the collapse to a complete 

graph takes place for A < Ac ^ N^''. In both cases (i.e 
for 6 > 3/2) for any A < 1, it is always possible to take N 
large enough to make the decay rate XN^ large enough 
so that finite size effects can be neglected. But when 
6 < 3/2 this is no longer true because {k) ^ N and the 
decay rate of completely connected nodes remains finite 
(7 = 0) even when N ^ 00. Beyond a finite dissipation 
rate Ac, the network collapses to the complete graph. 
Fig. 13 fully confirms the theoretical insight discussed 
above. When A > Ac, where Ac ^ N ' (see inset), the 
dynamics reaches a network with finite average degree, 
whereas for A < Ac a collapse to the complete graph is 
observed, with a transition which is sudden and discon- 
tinuous. In the case in which the Uk = uwe also observe 
a transition from a finite average connectivity to a aver- 
age connectivity of order N. The transition occurs for 
values of Ac ~ but the transition is rather smooth. 
Having discussed the stationary state, let us go back to 
the Eq. ||SJ) for the avalanche size distribution. We focus 
on the region 1 ^ A ^ iV~^ where the network is not 
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FIG. 3: Data collapse of the avalanche distribution P{s) 
with s > 7 for networks oi N = 10^ nodes averaged over 
200 A'^ time steps after equilibration. Here b = 2.5 and 
A = 0.01,0.02,0.04,0.06,0.18,0.31. The data collapse was 
done taking the theoretical values r = 3/2 and ct = 2 of the 
exponents. Inset: the scaling function h{x) for the same data 
set. 

densely connected and dissipation effects are weak. An- 
ticipating that the avalanche size distribution acquires a 
cutoff Sc ^ A^'^, for some exponent a > 0, we postulate 
the scaling form P{s) ~ s~'^$(sA'^) with finite <I>(0) and 
^(x) — > faster than any power as x ^ oo. Such a scal- 
ing hypothesis is fully corroborated byjiumerical results, 
as shown in Fig. |31 This corresponds 15) to an analogous 
scaling form 

x(z)^i-(i-zr'i/.(^i^) (11) 

for the generating function for A <C 1 and 1 — z ^ y . 
Setting for convenience 1~ z — xX^ , asymptotic analysis 
for A ^ 1 shows that the leading orders of Eq. are 

(12) 

where 

/? = 2, forfe>2 (13) 

/3 = ^ for&<2. (14) 



Note that c ~ 1/|&— 2| diverges when 6 ^ 2 is approached 
from both sides. Dividing Eq. (|12|l by A°" and taking the 
scaling limit A — > with x finite, we find a non-trivial 
result with all three terms finite if wc choose 



T = 3/2, a = 2 for > 2 (15) 
T = 1 + 1/6, cr = b/{b - 1) for b<2. (16) 



and the scaling function is the inverse of x{h) = ^^/[l — 
c/i^]'^. In particular h{x) — > c^^^ for x ^ oo and 
h{x) ^ x^^^ for X <^ 1. For b > 2 the solution coin- 
cides with that of other mean field models 
[y^ c/x + 4— y/ c/x]/{2c) and perfectly matches numerical 
simulations for a range of values of A (see the inset of Fig. 
I^J. It is easy to check that the model with Uk decaying 
slower than 1/fc falls in the (3 = 2, t = 3/2 a = 2 uni- 
versality class. In conclusion we have shown how a slow 
growth dynamics and a fast relaxation through avalanche 
events can generate a dynamical network with given de- 
gree distribution. The stationary state is critical in the 
sense that the avalanches of all sizes occur, and is reached 
spontaneously without fine tuning of external parame- 
ters as long as the dissipation rate is larger than a given 
threshold Ac ^ 7V~^. For smaller dissipation rates, the 
network collapses to the complete graph. While the de- 
tailed solution depends on the particular simplicity of 
the model chosen, the generic picture may apply to a 
wider class of systems and capture some features of the 
non-linear and intermittent behavior of real systems. We 
acknowledge F. Vega-Redondo for useful comments and 
discussions. This work is supported in part by the Euro- 
pean Community's Human Potential Programme under 
contract HPRN-CT-2002-00319, STIPCO. 
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